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ABSTRACT 

A commonly  occurring  class  of  integral  equations  which  arise  regularly  in 
applications  are  the  separable  first  kind  Abel-type  integral  equations  such  as 

CD 

(*)  s(y)  = / k^  (y)  (x)  (x2  - y2)  Wu(x)dx,  0 < u < 1,  x 1 Y i 0 • 

y 

Though  the  mathematical  properties  of  this  equation  (such  as  conditions  for  the 

existence,  uniqueness  and  smoothness  of  its  solutions,  its  improperly  posed 

nature,  the  existence  of  inversion  formulas,  etc.)  have  been  examined  in  some 

detail  in  the  literature,  its  numerical  solution  poses  a number  of  difficulties 

especially  when  s(y)  is  only  available  as  discrete  observational  data 

(d.)  = {d.  = s(x.)  + e.,  i = 1,2, ...,n;  c.  discrete  random  errors)  . 
1111  l 

In  this  paper,  we  (i)  review  the  numerous  applications  in  which  the  solutions  of 
equations  like  (*)  for  discrete  observational  data  (d.)  is  the  basic  step, 

(ii)  with  respect  to  given  discrete  observational  data  { d ^ } , compare  the  use 
of  pseudo-analytic  methods  and  the  direct  evaluation  of  its  inversion  formulas 
as  a basis  for  solving  (*) , (iii)  propose  a specific  algorithm  based  on  the 
conclusions  of  (ii),  and  (iv)  examine  the  consequences  of  the  fact  that,  for  <*), 
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linear  functionals  defined  on  its  solution  u(x)  can  be  redefined  as  linear 


functionals  on  the  data  s(y).  The  justification  for  the  latter  is  that,  in 


applications  involving  separable  first  kind  Abel-type  integral  equations 


inferences  are  usually  based  on  (linear)  functionals  defined  on  u(x) 


u(x)  itself.  This  point  is  illustrated  with  an  example  from  metallurgy 
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significance:  and  explanation 


In  qeneral,  it  is  not  possible  to  determine  three-dimensional  structure 
(such  as  the  distribution  of  carbon  particles  in  steel,  the  size,  shape  and 
position  of  weld  defects  in  welded  thick  steel  plates,  and  the  distribution  of 
cells  in  and  structure  of  biological  and  zoological  specimens)  directly.  It  is 
necessary  to  make  indirect  observations  such  as  are  obtained  by  takinq  X-ray 
photographs  (two-dimensional  projections)  of  the  three-dimensional  structure,  or 
by  cutting  open  (taking  (random)  two-dimensional  (planar)  sections  of)  the 
three-dimensional  structure  and  photoqraphing , or  by  tnakinq  microscope  slides 
((small)  thin  sections  of  the  three-dimensional  structure)  for  examination  in  a 
(electron)  microscope.  Thus,  it  is  necessary  to  determine  information  (properties) 
about  the  original  three-dimensional  structure  from  the  two-dimensional  informa- 
tion (data)  obtained  from  the  indirect  observations  (such  as  the  X-ray  photo- 
graphs, photographs  of  planar  sections,  or  (electron)  microscope  photoqraphs) . 

In  a number  of  situations  (for  example,  when  the  three-dimensional  structure 
involves  axial  or  spherical  symmetry),  the  mathematical  formalism  which  relates 
properties  of  the  three-dimensional  structure  u(x)  to  corresponding  two- 
dimensional  data  s(y)  is  the  following  separable  first  kind  Abel-type  integral 
equation 

no 

(*)  s(y)  = / k (y) k2 (x)  (x2  - y2)  gu(x)dx,  0 < p < 1,  x>y^O, 

Y 

where  k^ (y) , k^(x)  and  g are  known.  For  example,  if  the  three-dimensional 
structure  consists  of  randomly  dispersed  (approximately)  spherical  carbon 
particles  in  steel,  and  u(x)  is  the  size  (radius  or  diameter)  distribution  of 

nr 

these  spt  al  particles,  then,  with  k (y)  - y/m,  m - / xu(x)dx,  k ( x ) = 1, 

1 0 

and  g = ) corresponds  to  the  random  spheres  model  of  stereoloqy  with 

s(y)  the  size  distribution  of  the  circular  sections  of  the  spherical  carbon 


p 

I particles  on  random  plane  sections.  In  fact,  the  mathematical  formalisms  for 

i a large  range  of  applications  are  defined  by  (*)  . A comprehensive  classifica- 

tion for  them  is  given  in  the  paper. 

Thus,  in  the  above-mentioned  situations,  the  determination  of  the  required 
three-dimensional  property  u(x)  reduces  to  solving  (*)  for  some  given  s(y). 
But,  in  application,  s(y)  is  not  usually  known  exactly,  but  is  only  aval?  ble 
as  discrete  observational  data 

d^  = s(x.)  + f , i = 1,2, ...,n  , 

where  the  c.  denote  discrete  random  errors.  The  paper  therefore  examine  in 
some  detail  the  question  of  how  to  solve  (*)  in  such  situations.  In  addition, 
because  of  its  practical  numerical  consequences , the  evaluation  of  linear 
properties  (functionals)  (such  as  mean  surface  area,  mean  volume,  etc. ),  defined 
on  the  solution  u(x) , of  the  form 

GO 

/ 6 (x)u  (x)dx  , 0 <_  y < » , 

y 

where  9 (x)  is  known,  is  discussed  in  some  detail  with  respect  to  a specific 
example  from  metallurgy  (namely,  the  determination  of  the  impact  strength  of  a 
given  consignment  of  steel) . 
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PREFACE 


This  paper  is  the  basis  for  the  invited  talk  given  at  the  Seminar 
organized  by  the  North-West  Division  of  SIAM  to  honour  the  retirement  of 
Professor  Arvid  T.  Lonseth,  Mathematics  Department,  Oregon  State  University, 
Corvallis.  It  was  held  in  Seattle  on  August  15  and  16  in  parallel  with  the 
81st  Summer  Meeting  of  the  American  Mathematics  Society. 

Support  from  the  Computer  Centre,  ANU  (Australian  National  University), 
and  the  Mathematics  Research  Center,  University  of  Wisconsin-Mad ison,  which 
ensured  that  participation  would  be  possible,  is  gratefully  acknowledged. 

In  compiling  a review  of  Abel-type  integral  equations  in  applications, 
it  has  been  necessary  to  discuss  various  points  with  colleagues  from  different 
disciplines  to  whom  the  author  is  indebted.  They  include  Drs.  Frank  de  Hoog 
(Computer  Centre),  John  Cleary  (Earth  Sciences),  Daryl  Daley  (Statistics), 

Pam  Davey  (Statistics),  Roger  Miles  (Statistics)  all  of  the  ANU  and 
Professor  J.  E.  Hilliard  (Material  Sciences)  at  Northwestern  University.  The 
author  wishes  to  record  particular  thanks  to  Tony  Jakeman,  CRES,  ANU.  The 
writing  of  this  report  has  benefited  substantially  from  lengthy  discussions 
with  him  about  numerous  finer  points  and  detail. 


Si.  Introduction:  What  is  an  Abel-type  Integral  Equation? 


Integral  equations  of  the  fora 
Y 

s(y)  ■ | K(y,x)  u(x)dx,  0 t y ( x i Y < -,  (1.1) 

y 

(and  its  equivalent 

f9 

s(y)  - j K(y.x)  u(x)dx,  0 < x < y < Y e - 
0 

with  y ■ Y - y,  x * Y - x)  are  known  as  /tree  kind  Volterra  integral  equation* 
with  kernel  K(y,x).  The  kernel  K(y,x)  is  said  to  be  non-singular  or  singular 
depending  on  whether  K(y,x)  is  continuous  or  discontinuous  on  the  triangular 
region  0 ( y < x ( Y.  In  the  special  situation  where  K(y,x)  takes  the  f"rw 

K(y.x)  - k(y,x)  / (xp  - yP)p.  0 < u < I.  p > 0.  (1.2) 

with  k(y,x)  continuous  on  0(y(x(Y  and  k(y,y)  t 0,  the  equations  (1.1) 
are  called  generalized  first  kind  Abel  integral  equations  or  first  kind  Abel-  ,-pe 
integral  equations. 

The  reason  for  the  name  is  easily  explained.  When  he  posed  and  solved  in 
1823  the  falling-particle  problem  (see  Lonseth(1977))  as  a generalisation  of 
Huygen's  problem  of  the  isochronous  pendulum  (1673),  Abel  obtained  the  equation 

| u(x)  (y  - x)'*1  dx  « T(y)  (1.3) 

0 

which  now  bears  his  name,  where  T(y)  denotes  the  time  for  a particle  to  fall 
without  friction  a distance  y,  g the  Earth's  gravitational  constant,  and  u(x) 
the  path  the  particle  must  be  constrained  to  fall  in,  if  its  time  of  fall  is  to 
coincide  with  some  prescribed  function  T = T(y)  of  the  distance  fallen.  The 
structural  similarity  between  (1.1),  with  K(y,x)  defined  by  (1.2).  and  (1.3)  is 
obvious.  In  fact,  the  similarity  goes  deeper  than  this  (see,  for  example, 

Atkinson  (1974))  and  is  the  main  reason  for  the  definitions.  The  pertinence  of 
this  point  will  become  clearer  in  the  sequel. 

As  we  shall  see  below,  a very  important  subclass  of  first  kind  Abel-type 
integral  equations  is  generated  when  k(y,x)  is  separable,  namely 

k(y,x)  * k, (y)  k2(x).  (1.4) 

For  example,  when  (1.1)  takes  the  form 

00 

s(y)  * | k,(y)k2(x)  (x  - y)~w  u(x)dx,  0<p<l,  x>y>0,  (1.5) 

y 
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u(x)  is  defined  by  the  following  inversion  formulas,  when  u(x)  k,  (x)  is 
continuous.  k2(x)  t 0,  x e [0,-),  and  s(y)/k,(y)  is  differentiable. 


- - Iwka(x)]'1  sin  (px)  i { ( (y  - x)v'1  dyj  . 

x 


(1.6a) 


(l .6b) 


We  shall  refer  to  this  important  subclass  as  separable  /tret  kind  Abel-type 
integral  equations. 

Since  the  addition  of  the  term  Xu(y),  X»*0,  to  the  right  hand  side  of 
(1.1)  converts  it  to  a second  kind  Volterra  equation,  namely 

Y 

s(y)  ■ Xu(y)  ♦ | K(y,x)u(x)dx,  0 < y < x < Y < «,  (1.7) 

y 

we  shall  refer  to  the  corresponding  equations  as  second  kind  Abel-type  integral 
equations,  when  the  kernel  K(y,x)  in  (1.7)  takes  the  form  (1.2),  and  separable 
second  kind  Abel- type  integral  equations  when  the  kernel  K(y,x)  in  (1.7)  takes 
the  form  (1.2)  and  k(y,x)  is  separable. 

Though  the  Abel -type  integral  equations  which  arise  most  frequently  in 
applications  are  of  the  first  kind,  the  second  kind  are  also  found.  In  particular, 
they  arise  in  the  analysis  of  thin  sections  in  transmission  electron  microscopy 
(e.g.  Goldsmith  (1967)).  The  most  commonly  occurring  Abel-type  integral  equations 
in  applications  are  certainly  the  separable  first  kind  equations. 


2 


f 


11.  Introduction:  What  is  an  Abel-type  Integral  Equation? 


Integral  equations  of  the  form 
Y 

*(y)  ■ | K(y.x)  u(x)dx,  0 ( y ( x ( Y < ».  (1.1) 

y 

(and  its  equivalent 

*(y)  “ | K(y,x)  u(x)dx,  0«x<y<Y<-, 

0 

with  y = Y - y,  x « Y - x)  are  known  as  first  kind  Vcl  terra  integral  equation* 
with  kernel  K(y,x).  The  kernel  K(y,x)  is  said  to  be  non-singular  or  singular 
depending  on  whether  K(y,x)  is  continuous  or  discontinuous  on  the  triangular 
region  0 < y i x s Y.  In  the  special  situation  where  K(y,x)  takes  the  form 

K(y»x)  ■ k(y.x)  / (xp  - yp)M,  0 < u < 1,  p > 0,  (1.2) 

with  k(y,x)  continuous  on  0 ( y S x £ Y and  k(y,y)  t 0,  the  equations  (1.1) 
are  called  generalized  first  kind  Abel  integral  equations  or  first  kind  Abel- type 
integral  equations. 

The  reason  for  the  name  is  easily  explained.  When  he  posed  and  solved  in 
1823  the  falling-particle  problem  (sec  Lonseth(I977))  as  a generalization  of 
Huygen’s  problem  of  the  isochronous  pendulum  (1673),  Abel  obtained  the  equation 

J u(x)  (y  - x)'*1  dx  « /2f  T(y)  (1.3) 

0 

which  now  bears  his  name,  where  T(y)  denotes  the  time  for  a particle  to  fall 
without  friction  a distance  y,  g the  Earth's  gravitational  constant,  and  u(x) 
the  path  the  particle  must  be  constrained  to  fall  in,  if  its  time  of  fall  is  to 
coincide  with  some  prescribed  function  T = T(y)  of  the  distance  fallen.  The 
structural  similarity  between  (1.1),  with  K(y,x)  defined  by  (1.2).  and  (1.3)  is 
obvious.  In  fact,  the  similarity  goes  deeper  than  this  (see,  for  example, 

Atkinson  (1974))  and  is  the  main  reason  for  the  definitions.  The  pertinence  of 
this  point  will  become  clearer  in  the  sequel. 

As  we  shall  see  below,  a very  important  subclass  of  first  kind  Abel -type 
integral  equations  is  generated  when  k(y,x)  is  separable,  namely 

k(y,x)  = k, (y)  k2(x).  (1.4) 

For  example,  when  (1.1)  takes  the  form 

00 

s(y)  = | k,(y)k2(x)  (x  - y)*p  u(x)dx,  0<p<l,  x>y>0,  (1.5) 

y 
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u(x)  is  defined  by  the  following  inversion  formulas,  when  u(x)  k,  (x)  is 
continuous,  kj(x)  / 0,  x e I0,“>),  and  s(y)/k,(y)  is  differentiable. 


u(x)  - - Ink  (x)  J~ * sin  u*  ^ — — 

LkJyHy-x)1  U I y .. 

•r— 

X 


- - lTTk,(x)]‘‘  sin  (mO^j  j (y  - x)W_l  dyj  . 

x 

Me  shall  refer  to  this  important  subclass  as  separable  first  kind  Abel-type 
integral  equations. 

Since  the  addition  of  the  term  Au(y),  Xyo,  to  the  right  hand  side  of 
(1.1)  converts  it  to  a second  kind  Volterra  equation , namely 

Y 

s(y)  » Xu(y)  ♦ | K(y,x)u(x)dx,  0 < y < x < Y < -,  (1.7) 

y 

we  shall  refer  to  the  corresponding  equations  as  second  kind  Abel-type  integral 
equations,  when  the  kernel  K(y,x)  in  (1.7)  takes  the  form  (1.2),  and  separable 
second  kind  Abel- type  integral  equations  when  the  kernel  K(y,x)  in  (1.7)  takes 
the  form  (1.2)  and  k(y,x)  is  separable. 

Though  the  Abel-type  integral  equations  which  arise  most  frequently  in 
applications  are  of  the  first  kind,  the  second  kind  are  also  found.  In  particular, 
they  arise  in  the  analysis  of  thin  sections  in  transmission  electron  microscopy 
(e.g.  Goldsmith  (1967)).  The  most  commonly  occurring  Abel-type  integral  equations 
in  applications  are  certainly  the  separable  first  kind  equations. 


(16a) 


(1.6b) 


5 2.  Preliminaries  and  Notation 


Because  it  is  impossible  to  cover  in  this  paper  all  the  major  aspects 
associated  with  the  theory  and  application  of  Abel-type  equations,  we  cite 
briefly  source  information  for  material  not  explicitly  discussed  in,  but  closely 
related  to  the  context  of,  the  paper,  as  well  as  define  notation  anJ  concepts 
required  subsequently. 

2.1  Existence  and  Uniqueness  of  Solutions  of  Abel-type  integral  Equations. 

Atkinson  (1974)  has  shown  that  an  existence  theory  can  be  constructed 
for  Abel-type  equations  (in  the  form  equivalent  to  (1.1))  when  the  kernel  K(y,x) 
has  the  very  general  form 

K(y.x)  = k(y,x)  / (yP  - xP) , 0 < x < 1,  p > 0, 

with  k(y,x)  continuous  on  0 ( x ( y ( Y and  k(y,y)  f 0.  In  addition,  the 
paper  is  a source  reference  to  previous  work  on  the  subject. 

2.2  The  Probabilistic  Basis  of  Stereology 

Excellent  discussions  of  the  geometric  probability  basis  of  the  theory 
of  estimating  the  properties  of  the  distribution  of  one  material  in  another  by 
investigating  the  observed  patterns  on  plane  and  line  intersections  (stereology) 
are  now  available: 

(i)  Moran  (1972)  for  a lucid  introduction. 

(ii)  Miles  and  Davy  (1977a),  (1977b),  Davy  and  Miles  (1977)  and  Davy  (1977) 

for  a rigorous  examination  of  the  probabilistic  assumptions  and  structure. 

(iii)  Santalo  (1953),  (1955)  for  a lucid  introduction  to  the  mathematical 

foundations  of  geometric  probability  (namely,  integral  geometry)  and  for 
one  of  the  key  examinations  for  convex  particles. 

2.3  Preliminaries  and  Notation  for  the  Estimation  of  Size  Distributions 


For  some  given  three-dimensional  object  (such  as  a sample  of  steel, 
cancer  infected  tissue,  a conglomerate  of  lipid  cells,  etc.),  the  aim  is  to  derive 
information  about  the  object's  three-dimensional  structure  from  lower  dimensional 
observations  of  it  (such  as  is  obtained  with  linear  probes  as  well  as  with  planar 
and  thin  sections).  For  example,  if  the  phase  of  interest  in  the  object  is 
particulate  (such  as  lipid  cells  in  a conglomerate,  cancer  cells  in  tissue  or 
inclusions  in  steel),  deductions  about  structure  for  application  purposes  may 
involve  estimation  of  some  or  all  of  the  following  properties  of  that  particulate 
phase:  shape,  number,  volume,  size,  surface  area  and  length.  When  simple  logic 
is  unable  to  yield  the  required  information  for  subsequent  inferences,  mathematical 
methods  based  on  geometric  probability,  statistics  and  numerical  analysis  are 
required  to  formulate  and  solve  the  equations  relevant  to  assumptions  made  about 
the  global  geometry  of  the  object. 

In  the  sequel,  because  our  interest  is  in  Abel-type  integral  equations,  we 
restrict  attention  to  objects  consisting  of  two-phases,  a and  6,  with  the 
ct-phase  made  up  of  similarly-shaped  convex  particles  and  the  6-phase  the  matrix 
in  which  the  a-phase  is  embedded.  It  is  assumed  that  the  a-phase  particles  are 

(i)  convex, 

(ii)  similarly  shaped  (i.e.  definable  in  terms  of  only  one  parameter  x) , and 

(iii)  randomly  distributed  in  a convex  field  (i.e.;  the  shape  of  the  object 
itself  is  also  convex). 


I 


We  denote  these  particles  by  {Px},  with  parameter  x,  and  assume  that  their 
probability  site  density  is  g(x)  * G'(x)  (G(-«>)  * 0). 

In  addition,  let  the  probability  density  of  the  observations  be 
(a)  ■ Z'(a)  where  a denotes  either  area  or  length,  depending  on  particle 
shape  and  the  type  of  observation  process  implemented.  Write  an  as  the  maximum 
over  the  sizes  a.  Also  dependent  on  these  two  factors,  use  of  the  notation  V 
or  L will  be  required:  let  V be  the  volume  of  the  particle  P,  (that  is,  the 
convex  particle  with  size  parameter  x * 1),  S its  surface  area  and  C its 
integral  of  mean  curvature.  When  these  quantities  are  not  appropriate,  we  use  L 
to  denote  length  of  P,.  Denote  by  the  average  number  of  particles  per  unit 

length  of  linear  probe,  by  the  average  number  per  unit  area  of  plane,  and  by 

Ny  the  average  number  per  unit  volume  of  convex  field.  When  this  convex  field  is 
transparent,  we  require  its  volume  Vp  and  surface  area  Sp. 


Product  Integration 


In  developing  numerical  methods  for  the  direct  evaluation  of  inversion 
formulas  such  as  (1.6a)  and  (1.6b),  it  is  necessary  to  evaluate  integrals  of  the 
form 


1(F)  •=  F(y)dy  = w(y)f(y)dy  = I(w;f)  , 0 < x < X < » , 


where  the  notation  allows  for  the  possibility  that  the  integrand  F(y)  can  be 
factorized  into  a smooth  part  f(y)  and  a non-smooth  (ill-behaved)  part  w(y). 

In  the  application  of  a standard  numerical  quadrature  procedure  to  1(F),  the 
strategy  is  to  approximate  F(y)  by  some  member  of  a family  of  smooth  functions 
{Fn(y)}  (which  depend  on  some  parameter  n)  all  of  which  can  be  integrated 
analytically.  The  most  commonly  used  family  is  the  polynomials  of  degree  n 
which  interpolate  F(y)  at  n + 1 points  on  [x,X] , viz.  (pn(F;y)}. 

This  strategy  is  implemented  independently  of  whether  or  not  F(y) 
possesses  a factorization  into  a product  of  a smooth  and  non-smooth  component. 
Though  convergence  for  Riemann  integration  can  be  proved  for  integrals  which 
possess  such  a factorization,  as  long  as  the  non-smooth  component  is  itself 
integrable,  this  does  not  guarantee  the  convergence  of  more  sophisticated  quadra- 
ture rules  than  Riemann  integration.  See,  for  example,  Elliott  (1977)  and 
the  historical  references  mentioned  there.  In  addition,  even 
if  such  a method,  which  does  have  guaranteed  convergence  for  the  integral  of 
interest,  is  used,  the  efficiency  associated  with  its  implementation  may  be 
grossly  inferior  when  compared  with  its  use  to  evaluate  an  integral  with  a 
smooth  integrand. 

The  resulting  conclusion,  which  appears  to  have  been  initially  identified 
and  implemented  by  Young  (19S4),  is  that  the  strategy  for  the  construction  of 
numerical  quadrature  procedures  for  such  integrals  should  take  into  account 
explicitly  the  fact  that  such  factorizations  of  the  integrand  F exist.  The  aim 
is  then  to  choose  a factorization  of  F(y)  into  smooth  and  non- smooth  parts,  f(y) 
and  w(v),  respectively,  in  such  a way  that  the  chosen  family  of  smooth  functions 
{fji  (y)J,  from  which  an  approximation  to  f is  taken,  are  such  that  the  integrals 


I(w;fn(y)) 


' w(y)fn(y)dy 


can  be  evaluated  analytically,  thereby  integrating  out  the  non-smooth  component 
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w(x)  exactly.  The  more  common  choice  of  the  (fn(y)}  are  linear  combinations 
of  basis  (coordinate)  functions  (4k  (y)  } , viz. 

f(y)  - 5 a.(n)  <My)  , (2. 

i=l  1 1 


such  that  the  integrals 


A 

I(w;<j>i(y))  = J w(y)$.(y)dy 


can  be  evaluated  analytically.  Thus,  for  a given  F(y),  after  an  appropriate 
choice  of  w(y)  and  {<J>i(y)}  have  been  made  such  that  the  integrals  (2.4)  can 
be  evaluated  analytically,  the  problem  of  product  integration  is  reduced  to  an 
examination  of  appropriate  facets  of  approximation  theory  concerning  the 
approximation  of  some  smooth  f(x)  by  functions  fn(y)  of  the  form  (2.3). 

The  justification  for  introducing  product  integration  as  the  strategy 
for  evaluating  integrals  of  the  mentioned  type  can  be  based  not  only  on  its 
success  as  a strategy  for  solving  singular  (integrable)  integral  equations  (see, 
for  example  Weiss  and  Anderssen  (1972),  Weiss  (1972),  and  Anderssen,  de  Hoog  and 
Weiss  (1973))  as  well  as  evaluating  singular  (integrable)  integrals  (see,  for 
example,  de  Hoog  and  Weiss  (1973)  and  Atkinson  (1975),  Part  I,  Chapter  5),  but 
also  on  the  fact  that,  in  statistical  contexts  such  as  the  evaluation  of 
functionals  on  particle  size  distribution,  product  integration  estimators  some- 
times have  superior  properties  when  compared  with  the  standard  statistical 
estimators  (see,  for  example,  Anderssen  and  Jakeman  (1975a)). 


2.5  Spectral  Differentiation 


As  we  shall  see  below,  the  key  to  obtaining  stable  computational 
strategies  for  the  direct  evaluation  of  inversion  formulas  (such  as  (1.6a) 
and  1.6b))  on  discrete  data  is  a stable  procedure  for  the  numerical 
differentiation  of  observational  data,  which  can  be  formulated  as  follows: 
for  given  observational  data 


d.  = s(x.)  ♦ e.  , 0 = x < x < 
j y j ' o » 


- < xn  - 1. 


s(x)  = H(x  - y)  u(y)  dy  = Lu  , say  , 


where  H(z)  denotes  the  Heaviside  unit  step  function.  Because  numerical 
differentiation  is  improperly  posed  in  that  small  perturbations  in  the  data 
can  yield  large  perturbations  in  the  solution,  standard  procedures  based  on 
an  abstract  formalism  (such  as  finite  difference  methods)  are  inappropriate 
in  that  they  function  independently  of  the  signal/noise  structure  in  the 
observational  data  (2.5)  to  which  they  are  applied.  It  is  therefore  necessary 
to  appeal  to  sophisticated  concepts  and  tools  to  obtain  stable  procedures 
which  do  allow  a separation  of  the  signal  from  the  noise  in  observational  data 
to  be  differentiated. 


Two  possible  strategies  are: 


2.5.1  Regularization  {Cullum(197l) , Anderssen  and  Blooufield(1974b) 
and  Anderssen (1 977 ) } For  the  observational  data  (2.5),  estimate  the  derivative 

of  the  signal  s(x),  namely  S(x),  as 


u(a;xk),  k * 0,1,2,  ....  1C,  0 « x0  < x,  < x,  < — < x^  » 1,  (2.7) 

where  u(a;x)  is  the  function,  contained  in  the  Sobolev  space  for  which 

u { III"  - *ll!  ♦ ° IMIJ  }»  (2.8) 

when  evaluated  on  the  data  (2.5),  attains  its  minimum,  where  ||*||  denotes  the 
norm  in  W^',  namely  ** 


IN 


2 

q 


I lull?  - ||«<» 


2 

0 


with 

« d*u/dxA,  l * 1,  2,  ...»  q, 
and  j 

iiuii;  ■=  | u*  dx. 

0 


2.5.2  Wiener  Filtering  (Hannan  (1970),  Anderssen  and  Bloomfield  (1974a) 
(1974b),  Koopmans  (1974)  and  Bloomfield  (1976)}  Assuming  the  observational  data 
(2.5)  has  been  generated  by  stationary  stochastic  processes,  estimate  s(x)  as 


u(£(w)  ; x.)  = £ (iNw  ) exp  (ikw.)  I (w  ) v (w.),  1 

-1T<W.<TI  J J J J , 

3 


k = 0,1,2 K,  xk  - (k-l)A,  A * 1/K,  i - /T,  , 

A 

where  l (w)  is  the  chosen  approximation  for  the  spectral  window 


(2.9) 


U*0  - 8s(w)  / (gs(w)  ♦ ge  (w)) 


(2.10) 


with  gs(w)  and  ge(w)  denoting  the  spectral  density  functions  of  s(x)  and 
{Cj},  respectively,  where 

, K 

v(w.)  « (K  ♦ 1)  l d.  exp  (-  ikw.) 

} k=0  J 

defines  the  finite  Fourier  coefficients  of  the  data,  and  where  w^  « 2uj/K.  # 
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Computationally,  both  suffer  from  the  same  defect.  The  former  will 
only  yield  a satisfactory  approximation  when  an  appropriate  value  has  been 
estimated  for  a from  the  data,  while,  for  the  latter,  an  appropriate 
estimate  4(w)  of  S.(w)  is  required.  In  both  situations,  estimates  of  the 
relevant  unknown  can  be  evaluated,  if  some  appropriate  criterion,  which 
reflects  the  signal/noise  structure  in  the  given  data  (2.5),  is  used.  Kahba 
(1976),  (1977)  has  investigated,  in  a quite  general  context,  the  use  of  cross- 
validation  for  estimating  a,  while  the  statistical  and  time  series  literature 
contains  numerous  non-parametric  approximations  t(w)  for  i.(w)  (see,  for 
example,  Jenkins  and  Watts  (1968)  and  Bloomfield  (1976)}. 

On  the  other  hand,  Anderssen  and  Bloomfield  (1974a, b)  have  proposed  a 
parametric  model  for  2(w).  It  is  based  on  the  fact  that,  under  appropriate 
conditions  (including  d0  = d^  = 0)  and  stationarity  assumptions  about 
(dj),  if  I(w)  is  modelled  as 

q 

(w)  = [1  ♦ a(w/A)*  l (w/A)*j]  (2.11) 

j=0 

then 

u(a;xk)  = u(£(w);  xk)  , k = 1,  2,  ....  K.  (2.12) 

Not  only  does  this  equivalence  yield  a basis  for  explaining  why  the  regulariza- 
tion (2.8)  works,  it  also  yields  (by  assuming  that  it  is  something  special  and 
therefore  should  be  utilized),  in  conjunction  with  the  work  of  Whittle(19S2)  on 
estimating  parametric  models  for  gs(w)>  a basis  for  estimating  a with  a 
specific  optimal  sense.  In  fact,  if  it  is  assumed  that  the  errors  (ej)  are 
uncorrelatcd  with  mean  zero  and  variance  27ib  so  that  ge(w)  = b,  then  the 
parametric  model  for  gs(w)  becomes 


gs(w)  * b*W  (w)/(l  - X<q)  (w)). 


(2.13) 


Since  the  distribution  of  a stochastic  process  is  determined  by  its 
spectrum,  at  least  for  Gaussian  data,  a can  be  estimated  from  the  statistical 
properties  of  the  data  (dj).  In  fact,  for  fixed  q,  Anderssen  and  Bloomfield 
(1974a)  showed  that,  using  Whittle  (1952),  the  evaluation  of  a "maximum 
likelihood"  estimate  for  a reduces  to 


‘ a"  [~(K+l)/2  log{  l I(w  ) [1  ♦ {[  X <q)  (w)  | *-  1 }"*  l*1} 
L '•-iKWjSn  J J > 

* l log  [1  ♦ {[x<q)  (w)]  1 }•*  ] 1 

-TT<WjCTT  ( 1 -1 


(2.14) 


where  I(w^)  denotes  the  periodogram  defined  by 


1{wj)  = 2nW  | J0  dk  exP  (ikV 
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When  making  an  examination  of  the  occurrence  of  some  special  equation, 
such  as  Abel-type  integral  equations,  in  applications,  the  usual  practice  is  to 
list  various  disciplines  such  as  physics,  biology,  seismology,  medicine,  geometric 
probability  etc,  and  record  for  each  listed  discipline  t he  different  applications 
of  the  special  equation.  Such  a classification  for  Abel -type  integral  equations 
is  given  in  Table  1 {see  Appendix  - Table  1 }.  It  is  not  as  informative  mathe- 
matically as  it  first  looks.  It  clearly  indicates  that  Abel-type  integral 
equations  often  arise  in  applications,  but  gives  no  insight  into  why  this  happens. 
That  is,  it  gives  no  insight  into  why,  for  some  given  application,  it  is  an  Abel- 
type  equation  which  arises  as  the  mathematical  formalism  rather  than  some  other 
type  of  integral  equation  or  even  a differential  equation. 

On  the  other  hand,  the  various  applications  of  Abel-type  integral  equations 
could  be  classified  on  the  basis  of  the  underlying  mathematical  rationale  for  its 
occurrence.  This  is  the  motivation  for  the  classification  of  Table  2 {see  Appendix  - 
Table  2}.  What  now  becomes  clear,  and  therefore  is  the  justification  for  the 
classification  presented  in  Table  2,  are  the  reasons  why  Abel-type  integral  equations 
arise  and  will  continue  to  arise  so  often  in  applications.  The  principal  ones, 
namely  rationales  1,  4,  6 and  8,  are  either  directly  or  indirectly  related  to 
properties  of  axially  and/or  radially  symmetric  systems.  Thus,  the  rationales  1, 

4,  6 and  8 cannot  be  viewed  as  being  completely  independent.  In  fact,  there  is  a 
close  connection  between  1 and  6.  This  relates  to  the  way  in  which  a uniformly 
random  plane  which  sections  the  spherical  particle  Px  is  defined.  It  depends  on 
the  initial  choice  of  a random  direction  followed  by  the  choice  of  a random  point 
on  the  orthogonal  projection  of  Px  with  respect  to  the  chosen  direction  - see 
Watson  (1971),  §2,  and  Miles  and  Davy  (1977a).  The  geometric  properties  of  the 
sections  of  the  sphere  of  1 are  similar  mathematically  to  the  corresponding 
properties  of  the  paths  of  projected  rays  through  the  cylinders  and  spheres  of  6. 

Clearly,  the  rationales  2,  3,  5 and  7 are  completely  independent  of  axial 
and  radial  symmetric  considerations.  In  3,  it  is  the  formal  structure  of  the 
relationship  between  potential  and  kinetic  energy  which  determines  the  Abel -type 
equation  structure,  while  in  5 it  is  the  nature  of  the  kernel  in  Cauchy  integral 
equations.  In  7,  it  is  the  fact  that  the  variance  of  a Brownian  motion  process 
is  proportional  to  time  of  travel.  As  for  2,  the  explanation  is  directly  connected 
with  the  existence  of  inversion  formulas  for  separable  Abel-type  integral  equations. 
This  point  is  basic  to  much  of  the  discussion  in  the  sequel. 


Even  in  an  hour-long  talk,  time  does  not  permit  one  to  discuss  in  detail 
all  these  applications.  We  therefore  restrict  attention  to  the  specific  applica- 
tions which  arise  in  geometric  probability.  The  justification  for  this  is  three- 
fold. The  least  known  applications  of  Abel-type  integral  equations,  but  potentially 
some  of  the  most  important,  arise  in  geometric  probability.  Lucid  discussions  of 
applications  in  areas  other  than  geometric  probability  are  better  known  and  more 
readily  available  (for  example  Lonseth's  (1977)  discussion  of  sources  and  applica- 
tion of  integral  equations  and  the  references  listed  in  Tables  1 and  2).  The 
geometric  probability  framework  naturally  yields  applications  of  Abel-type  integral 
equations  with  kernels  of  the  form 

k,(y)k2(x)  / (x  - y)a,  x ) y i 0,  0 < o < 1,  a / !j. 

Even  though  this  last  fact  has  no  mathematical  significance  for  the  study 
of  the  various  generalizations  of  Abel's  original  equation  which  have  been  proposed, 
it  does  represent  positive  independent  motivation  for  studying  the  generalizations 
associated  with  the  above  mentioned  kernels  rather  than  some  other  generalization. 


We  therefore  turn  to  an  examination  of  the  (mathematical)  rationale  behind 
the  generation  of  Abel-type  equations  in  geometric  probability.  Within 
the  framework  of  geometric  probability,  it  is  that  aspect  of  this  subject  which  is 
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connected  with  (Underwood,  de  Wit  and  Moore  (1976),  p.v) 

"the  exploration  of  three-dimensional  space  (structure) , when 
only  two-  (or  one-)  dimensional  sections  through  solid  bodies 
or  their  projection  on  a surface  (or  a line)  are  available.", 

which  is  the  major  source  for  Abel-type  integral  equations  in  applications. 

This  is  clearly  illustrated  in  Table  2.  The  inter-disciplinary  activity 
associated  with  formulating  theory  and  implementing  (computational)  methods  to 
achieve  the  above  goal  (viz.,  extrapolation  from  discrete  (noisy)  observations 
in  lower  dimensions  (one  (lineal)  or  two  (areal))  to  structure  in  higher  (three)) 
is  now  referred  to  as  stereology.  In  fact,  there  is  an  international  society 
called,  not  surprisingly,  the  International  Stereological  Society  (ISS),  which 
aims  to  stimulate  communication  between  the  different  researchers  concerned  with" 
such  problems,  whether  their  interest  be  mathematical  or  non-mathematical , 
theoretical  or  pragmatic,  computational  or  qualitative,  experimental  or  instru- 
mental. Though  the  Foundation  Meeting  of  the  Society  was  only  held  in  1961,  the 
existence  of  the  Society  has  helped  pull  together  into  a limited  number  of 
sources  basic  information  about  stereology,  its  theoretical  and  experimental 
foundations,  its  applications,  and  its  instrumentation.  For  those  interested  in 
pursuing  this  point  further,  the  following  sources  are  suggested: 


1.  W.L.  Nicholson  (Editor),  Proceedings  of  the  Symposium  on 
Statistical  and  Probabilistic  Problems  in  Metallurgy,  Seattle, 
Washington,  4-6  August,  1971,  in  Special  Supplement  to  Adv. 

Appl.  Prob.,  December,  1972. 

A series  of  articles  which  illustrate  the  theoretical  and 
experimental  foundations,  the  applications  and  the  instrument- 
ation of  stereology. 

2.  E.E.  Underwood,  R.  de  Wit  and  G.A.  Moore  (Editors)  Fourth 
International  Congress  for  Stereology,  National  Bureau  of 
Standards  Special  Publication  *431,  U.S.  Government  Printing 
Office,  Washington,  Jan,  1976. 

Proceedings  of  the  Fourth  International  meeting  of  the  Inter- 
national Society  for  Stereology  which  was  held  at  the  National 
Bureau  of  Standards,  Gaithersburg,  Maryland,  U.S. A.  on  September 
4-9,  1975. 

A number  of  invited  as  well  as  contributed  papers  on  all  aspects 
of  stereology. 

3.  E.E.  Underwood,  Quantitative  Stereology,  Addison-Wesley,  Mass., 
1970. 


A text  book  discussion  of  the  subject  along  with  a quite  compre- 
hensive set  of  references  to  applications. 


,4.  R.T.  De  Hoff  and  R.N.  Rhines,  Quantitative  Microscopy,  McGraw- 
Hill,  New  York,  1968. 

Contains  a scries  of  invited  articles  by  various  experts  in  the 
field.  A collection  of  articles  on  some  of  the  theoretical  and 
experimental  foundations  as  well  as  the  applications  and  instru- 
mentation. 


5.  Journal  of  Microscopy  which  is  the  official  Journal  of  the 
International  Stereological  Society. 


-10- 


6.  Publications  of  Proceedings  of  the  meetings  of  the  Inter- 
national Stereological  Society  previous  to  that  cited  Tn  2. 
above.  Further  details  from  Prof.  Anna-Mary  Carpenter, 
University  of  Minnesota,  Department  of  Anatomy,  Minneapolis, 
Minn.  55455,  U.S.A. 

7.  Other  references  cited  in  this  paper  such  as  the  papers 
mentioned  in  §2.2. 


Specifically,  we  examine  the  problem  of  estimating  the  size  distribution 
of  three-dimensional  convex  particles  {Px}  in  an  opaque  medium  from  lower 
dimensional  observations  of  them.  The  formalisms  generated  essentially  fall 
into  three  categories  related  directly  to  the  nature  of  the  lower-dimensional 
data  (Santalo  (1955)  and  Jakeman  and  Anderssen  (1975),  §2  and  Table  1)).  The 
three  categories,  as  defined  by  the  lower  dimensional  data,  are: 

(a)  intersections  of  convex  particles  by  random  plane 
sections  taken  through  an  opaque  field; 

(b)  intersections  of  flat  (no  thickness)  convex  particles 
by  random  plane  sections  taken  through  an  opaque  field; 

(c)  intersections  of  convex  particles  with  random  linear 
probes  taken  through  an  opaque  field. 

The  key  is  this.  If  we  let  K(a,x)da  be  the  probability  that  a particle 
P_  , when  sectioned  or  probed  randomly,  has  an  area  or  length  of  intersection 
which  takes  a value  between  a and  a+da  , then  g(x)  and  z(a)  are  related 
by  the  first  kind  Volterra  integral  equation 

oo 

I xpK(a,x)g(x)dx  = c.Zia),  (3.1) 

f(a) 

where  the  parameters  f(a),  p and  c are  defined  for  the  three  categories  (a), 
(b)  and  (c)  above  in  Table  3.  For  example,  for  convex  particles,  Santal6  (1955) 
obtains  (3.1)  by  equating  two  different  expressions  tor  the  average  number  of 
particle  areas  in  the  range  (a, a+da),  which  lie  on  random  sectioning  planes.  One 
is  defined  in  terms  of  K(a,x)  while  the  other  is  defined  in  terms  of  z(a). 

TABLE  3 


Table  for  the  Parameters  f(a),  p and  c of  Equation  (3.1) 


(the  constants 

c»  an>  Nv, 

etc  are  defined 

in  §2  ) 

Data  Category 

££■1 

JE_ 

c 

(a) 

(*/*«)  ** 

1 

2*Na/(CNv) 

(b) 

a/an 

1 

4Na/(LNv) 

(c) 

«/an 

2 

4NL/(SNy) 
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Thus,  the  problem  of  generating  the  various  formalisms,  reduces  to 
determining  explicitly  the  structure  of  K(a,x).  In  Table  4 {see  Appendix  - 
Table  4}t  various  forms  taken  ty  K(a,x)  and  the  resulting  forms  for  (3.1) 
are  listed  along  with  the  associated  size  distribution  estimation  problem. 

Even  though  the  two  applications  associated  with  (c)  correspond  to  (numerical) 
differentiation,  they  are  included  because  (numerical)  differentiation 
corresponds  to  the  degenerate  situation  where  u * 0 and  k(y,x)  separable 
in  (1.2). 

For  example,  consider  the  size  distribution  estimation  problem  for 
circular  sections  observed  on  random  plane  sections  taken  through  an  opaque 
field  of  spheres.  From  the  point  of  view  of  estimating  K(a,x)  da,  the  key 
property  of  similarly  shaped  convex  particles  is  that  "an  area  a on  a 
planar  section  through  the  particle  Pt  corresponds  to  an  area  x2a  on  a 
planar  section  through  the  particle  Px".  Consequently,  we  obtain 


K(a,  l)da  = K(ax2,  x)  d (x2a) 

and  hence,  on  writing  a = ax2, 

K(a,  x)da  = K(a/x2,  1)  da/x2.  (3.2) 

However,  for  a sphere  of  radius  R,  since  "the  probability  the  sectioned 
circle  has  a radius  in  the  range  r to  r + dr"  equals  "the  probability 
the  sphere  is  cut  at  a distance  from  its  centre  in  the  range  x to  x ♦ dx", 
which  in  turn  equals  "xdx  (uniformly  distributed)",  it  follows  that 

K(a,  1)  = (4ir(irR2  - a))  = (4*(an  - a))'** 

which  yields 

K(a,  p)  = (4irp2(anp2  - a))"**  , p = x/R,  , (3.3) 


with  the  largest  observed  particle  scaled  to  correspond  to  p = 1. 

Thus,  for  spheres,  (3.1)  becomes,  on  using  (3.3)  and  the  appropriate 
results  of  Table  3, 


g(x)dx 


'(a/a^  2(1r(anx2-a))Js 


2ttN. 

A 

CN„ 


which  in  turn  reduces  to 


g(x)dx  = « 


y/yn  " 


( V-y2)1*  X 


r h(y) 


(3.4) 
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z(a)da  * h(y)dy  with  z(a)  « h(y)/2ny,  where  h (y)  denotes  the  probability 
density  for  the  radii  of  the  spheres.  On  introducing  the  substitution  u * xyn, 
the  first  kind  Abel-type  equation  known  as  the  random  spheres  model  is 
obtained : 


54.  The  Numerical  Analysis  of  Abel-type  Integral  Equations 


The  foundations  of  the  numerical  analysis  of  any  special  class  of 
equations  (such  as  Abel-type  integral  equations)  rests  heavily  on  the  known 
results  of  mathematical  analysis  about  the  properties  of  such  equations.  Such 
properties  include  existence,  uniqueness  and  smoothness  conditions  for  their 
solutions,  conditions  under  which  closed  form  solutions  and  inversion  formulas 
can  be  constructed,  existence  of  different  (but  equivalent)  representations  for 
their  solutions,  existence  of  families  of  complete  functions  for  which  the 
inversion  of  subclasses  of  the  given  equations  can  be  constructed  analytically, 
and  the  properties  of  discretized  versions  of  the  given  equations.  For  the 
separable  Abel -type  integral  equations,  they  not  only  form  the  most  important 
subclass  of  Abel-type  from  the  point  of  view  of  the  regularity  with  which  they 
arise  in  applications,  but  also  in  terms  of  the  major  extent  to  which  properties 
of  the  above  type  are  known. 

From  the  point  of  view  of  the  present  discussion,  the  important  proper- 
ties of  separable  Abel-type  equations  are: 

(i)  Existence,  Uniqueness  and  SmoothneB8.  As  already  mentioned  in 
the  Preliminaries,  the  source  reference  is  Atkinson  (1974).# 

(ii)  Weakly  Improperly  Posedness.  Because  the  numerical  differentia- 
tion of  discrete  observational  data  is  improperly  posed  in  the  sense 
that  small  perturbations  in  the  data  can  yield  large  perturbations  in 
the  solution,  the  inversion  of  first-kind  Abel-type  equations  for 
discrete  observational  data  is  regarded  to  be  weakly  improperly  posed, 
since  it  corresponds  to  a fractional  differentiation  of  the  data. 
Perturbations  6u(x)  in  the  solution  of  the  Abel  equation  (T1.4)  of 
Table  1,  with  kj  (y)  = y/m  and  k2(x,y)  = 1,  with  m = /0a  xu(x)dx, 
corresponding  to  different  perturbations  6s (y)  in  the  data  are 
illustrated  in  Table  5.# 

(iii)  Existence  of  Inversion  Formulas.  Along  with  appropriate  assumptions, 
various  extensions  and  generalizations  of  the  correspondence  between  the 
Abel-type  equation  (1.5)  and  the  inversion  formulas  (1.6a)  and  (1.6b)  are 
known.  For  example,  when  (1.5)  is  replaced  by 

q ry  k2(x)u(x) 

y Pk,(y)  — r dx  = s(y)  , 0 < y < Y , pa  - B > 0 , (4.1) 

J0  (yp-  xK) 

with  kt(z)  i 0 and  k2 (z)  / 0 , 0 < z < Y , the  inversion  formula 


pa+6-1 

u(x)  = -1700"  ta  + x *(x)1 


with 

and 
♦ (*) 


p sinau 

(pa  + B) 

r -wi 

f ^ dz 

TT 

Lm°)J 

l 

p sinau  f 

2p+b 

r d 

f s(zx)  ] pa+B  [ 

s(zx) 

‘ * J 

0 

(l-zP)‘-a 

|_d(zx) 

l k,(zx)  J zx  1 

k,(zx) 

14- 

. dz 

k,(0)J  dZ  ' 

— Li, 


A 


TABU-:  5 

Exemplification  of  the  Weakly  Improperly  Posed 
Nature  of  First  Kind  Abel-type  Integral  Equations 


fts(y) 

6u(x) 

2y(a2-  x2)V*‘  , Re  v > 0 

Special  Cases 

V = 4 iry/m 

(s(y)  bounded  on  0 i y 4 a) 

2x(a2-  x1)'’* 

(u(x)  unbounded  on  0 ( x ( a) 

v = 1 ny(a2-  y2)*Vm 
(lim  s(y)  = 0) 

Y*  a 

2x 

(lim  u(x)  = 2a) 

Y*  a 

2.  ^sin1  (!i0tt) 

2x(a2-  X2)*12  (a(a2-  x2)) 

(t  = a2-  y2) 

(bounded  on  0 (,  y 4 a) 

(unbounded  on  0 i x i a) 

3-  r^v^j xtV  «)/■ 

2x(a2-  x2)v_1  exp(a(a2-  x2)).  Re  v > 0 

(t  = a2-  y2) 

can  be  derived  from  equations  (2.1),  (2.2)  and  (2.5)  of  Atkinson  (1974). 
In  fact,  as  long  as  conditions  are  satisfied  which  guarantee  existence, 
uniqueness  and  appropriate  smoothness  of  their  solutions,  all  first  kind 
Abel-type  integral  equations  with  separable  kernels  have  inversion 
formulas.# 


Existence  of  Complete  Systems  for  which  the  Inversion  is  Known 
Analytically . 


ft*  = s , ft  : S,  - Sj  , (4.2) 

Sj  and  Sj  are  given  functional  spaces,  denote  a densely  defined 
operator  representation  for  some  given  first  kind  Abel-type  integral 
equation,  and  let  {^(x)}  denote  a explicitly-known  coordinate  system 
which  spans  the  domain  of  A , J)  (A)  , such  that  the  system  {#.  = Ai^  } 
is  also  known  explicitly  and  spans  the  range  of  ft  , R (A)  , dense  in  S2 
For  such  systems,  if  “ 

s - l <t>.  , n < » , 

k»l 

then  the  inversion  can  be  performed  analytically  to  yield 

n . . 

U ■ l a£nJ  Ik  , n < « , 
k=l 
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(iv) 

Let 


where 


TABLE,  6 

Some  Invertible  Coordinate  Systems*  for  the  First  Kind  Abel-type 


Integral  Equation**  s(y)  = 2 


(x2-y2)2 


dx  , 0 ( y i x i a < «. 


Source 

Reference 


j (y)  = (a2-y2)^. 

v Minerbo  and 

♦fc(y)  - (a2-y2)k-‘.  k > 1.  \(x)  = (a2-x2)k'\  k > 1 ^ (1969>- 


f<0(x)  * 1/Tt 


- t6Tk(l-2t),  k > 0. 

2 Tk(z)  Chebyshcv  polynomial 
of  degree  k , t = a2-y2  , 

6 = const.  > - *s. 


V‘>  ■ * 


<jF2^‘kik,B+1;!i  * 


Piessens  and 

Verbacten 

(1973). 


.(y)  = ft. (a2-  y2) 


4-k(x)  = r(p)  Fk  (a2-  X2) 


where  f.  (x)  and  F.  (y)  corresponding  to  any  systems  of 

Erdelyi (1954) 

functions  in  a table  of  Rieman-Liouville  integral  transforms: 

/ 

Fk(y)  = [r(p)]"1j  fk(x)  (y  - X)p_l  dx. 


* Note  1.  The  corresponding  invertible  coordinate  systems  for  other  first 
kind  Abel-type  integral  equations  can  be  obtained  by  transformation. 

**  Note  2.  Similar  systems  exist  for  the  equation  s(y)  = £ J^°  — -q  dx. 
The  obvious  source  is  a table  of  Wcyl  integral  transforms  - Erdelyi  (1954). 


when  the  mapping  A : D (A)  ■*  R (A)  is  one-to-one.  Because  of  the 
indirect  connection  with  variational  and  projection  methods,  we  shall 
refer  to  the  dual  system  (ijjk  » ^k  = Aipk } as  an  A-  invertible 
coordinate  system. 

Table  6 contains  some  examples  of  A-invertible  coordinate  systems  for 
some  commonly  occurring  Abel-type  integral  equations.  It  is  pertinent  to 
note  that  such  explicit  systems  can  only  be  constructed  for  Abel-type 
equations  for  which  the  kernel  K(y,x)  of  (1.2)  is  known  explicitly.  The 
existence  of  such  systems  is  basic  to  the  use  of  pseudo-analytic  methods 
which  will  be  discussed  in  some  detail  below.# 
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In  part,  the  existence  of  this  extensive  information  base  about  the 
mathematical  properties  of  separable  Abel-type  integral  equations  explains  why 
the  literature  on  the  numerical  analysis  of  such  equations  is  not  only  extensive 
but  also  contains  a great  variety  of  mathematically  independent  proposals  for 
their  solution  numerically.  They  include  the  standard  finite  difference 
approaches  of  Linz  (1968),  Brunner  (1973),  (1974),  product  integration  procedures 
of  Young  (1954),  Linz  (1967),  Noble  (1970),  Weiss  (1972),  Weiss  and  Andcrssen  (1972), 
and  Anderssen,  de  lloog  and  Weiss  (1973),  the  pseudo-analytic  methods  of  Minerbo 
and  Levy  (1969)  and  Piessens  and  Verbaeten  (1973),  use  of  finite  difference 
methods  by  Nestor  and  Olsen  (1960),  to  evaluate  the  inversion  formulas  and  the 
proposal  by  Anderssen  (1973),  (1976)  to  use  the  spectral  differentiation  procedure 
of  Anderssen  and  Bloomfield  (1974a, b)  to  evaluate  the  inversion  formulas  and  its 
implementation  by  Anderssen  and  Jakeman  (1975a, b)  for  Abel-type  equations  in 
stereology.  Except  for  the  recent  paper  of  Baev  and  Glasko  (1976),  the  direct 
application  of  regularization  does  not  appear  to  have  been  examined. 

Most  of  the  numerical  methods  proposed  in  the  stereological  literature 
for  the  approximate  solution  of  Abel-type  equations  correspond  to  the  application 
of  (standard)  finite  difference  techniques  (such  as  ones  based  on  the  use  of 
rectangular  quadrature  (h'atson  (1971))  and  variants  of  product  mid-point  (Wicksell 
(1925),  Meisner  (1967),  Baudhuin  and  Berthet  (1967),  and  IVeibel  et  al  (1969))  in 
conjunction  with  a data  grouping  strategy  to  obtain  a discrete  approximation  to  the 
size  distribution  data  s(y)  from  given  observations.  Few  appear  to  have  known 
that  an  inversion  formula  existed  and  even  fewer  used  it  (Meisner  (1967),  Appendix, 
makes  indirect  use  of  it  when  deriving  equations  (A. 12)).  Some,  as  cited  above, 
appear  to  have  used  product  integration  techniques  without  explicitly  ackno’. ledging 
that  this  was  the  situation.  For  example,  the  formulas  on  p.101  of  Tallis  (1970) 
correspond  to  the  use  of  a product-trapezoidal  technique,  while  those  of  (13)  of 
Meisner  (1967)  correspond  to  product  midpoint.  Few  appear  to  have  explored  the 
possibility  of  working  directly  with  the  cumulative  distribution  generated  by  the 
observational  data.  A number  have  proposed  methods  of  a more  specialized  nature 
such  as  the  parameter  estimation  procedure  of  Keiding  et  al  (1972)  and  the 
successive  subtraction  algorithm  of  Saltikov  (1967);  but  substantive  justification 
for  using  them  does  not  appear  to  exist,  while  strong  reasons  for  not  using  them 
have  been  put  forward  by  Jakeman  (1975),  §3.3,  and  Anderssen  and  Jakeman  (1975). 

Another  key  reason  for  the  variety  of  such  proposals  is  the  existence  of 
applications  where  the  data  s(y)  is  known  analytically,  such  as  the  calculation 
of  boundary  crossing  probabilities  for  Brownian  motion  and  the  associated 
applications  mentioned  in  Table  1.  In  fact,  many  of  the  methods  for  Abel -type 
equations,  which  do  not  involve  any  form  of  stabilization,  can  be  applied  success- 
fully to  the  solution  of  this  subclass  of  problems,  if,  with  respect  to  the 
chosen  discretization  of  the  given  Abel-type  equation,  the  required  discretiza- 
tion of  s(y)  is  computed  "sufficiently"  accurately.  However,  space,  time  and 
context  does  not  permit  even  a cursory  discussion  of  all  aspects  of  the  numerical 
analysis  of  Abel -type  equations,  and  thus  the  exact  data  subclass,  in  particular. 

The  interested  reader  is  referred  to  Weiss(1972),  and  Anderssen,  de  Hoog,  and 
Weiss  (1973),  as  well  as  the  more  recent  work  of  Holyhead  and  McKee  (1976),  (1977). 

We  therefore  limit  attention  to  the  numerical  solution  of  first  kind 
Abel-type  equations  of  the  form 

Y 

s (y)  - | kj (y)  k?(x)  (x  - y)’Pu(x)  dx,  0 < p < 1,  0<y<x<Y<<*>,  (4.3) 

y 

where  p,  Y,  k^y)  and  k2(x)  are  known  explicitly,  and  when  the  data  s(y)  is 
only  available  as  discrete  observations: 

{d.}  = (d.  = sfx^  ♦ e^,  i = 1,  2,  ...»  n;  discrete  random  errors).  (4.4) 
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From  the  point  of  view  of  proposing  a numerical  procedure  for  the 
approximate  solution  of  (4.3)  for  the  given  data  (dj),  the  aim  is  to  obtain 
from  the  given  data  as  much  information  as  possible  by  maximizing  knowledge 
about  the  properties  of  the  solution  of  (4.3).  It  is  now  necessary  to  be 
sympathetic  with  the  "red  Queen's  philosophy" 

"I  could  have  done  it  in  a much  more  complicated 
way",  said  the  red  Queen,  imnensely  proud. 

Lewis  Carroll 

in  that,  in  trying  to  achieve  the  above  goal,  it  may  be  necessary  to  be  devious 
and  complicated  as  well  as  crafty. 

By  themselves,  parsimony  and  simplicity  should  not  now  be  the  guiding 
strategy.  Instead,  a suitable  trade-off  between  the  competing  claims  must  be 
made  and  implemented.  For  the  solution  of  first  kind  Abel-type  equations,  the 
competing  claims  are: 

(a)  smoothing  of  the  data  (separation  of  the  errors  (noise) 
from  the  structure  (signal)  in  the  data); 

(b)  construction  of  a stabilized  computational  strategy  for 
the  solution  of  the  formalism  (4.3). 

Many  of  the  methods  proposed  fail  to  do  either,  but  these  are  usually  the 
methods  proposed  for  exact  data,  whereas  the  above  claims  are  first  and  foremost 
a direct  consequence  of  the  assumption  that  the  data  is  discrete  and  observational. 
Additional  motivation  for  (b)  is  the  previously  mentioned  weakly  improperly  posed 
nature  of  first  kind  Abel-type  equations. 

This  weakly  improperly  posed  nature  is  also  the  reason  for  the  failure 
of  standard  finite  difference  methods  for  the  solution  of  (4.3)  when  the  data  is 
discrete  and  observational  (or  discrete,  but  not  computed  with  sufficient  accuracy). 
As  shown  in  Anderssen  (1976),  §2,  they  have  no  natural  stabilization  against  the 
propagation  of  small  errors  in  the  data  into  unacceptably  large  errors  in  the 
solution.  In  fact,  finite  difference  methods  fail  to  meet  either  of  the  claims 
(a)  or  (b) . Insight  into  why  this  is  so  can  be  gained  from  an  examination  of  the 
use  of  finite  difference  methods  for  the  differentiation  of  discrete  observational 
data,  since  the  solution  of  Abel-type  equations  corresponds  to  a fractional 
differentiation  of  the  data. 

There  are  two  obvious  ways  by  which  the  above  competing  claims  can  be 
satisfied.  They  are: 

1.  Apply  either  the  inversion  formula  (1.6a)  or  (1.6b)  to  filtered  data. 

The  simplest  implementations  of  this  strategy  are  the  pseudo-analytic  methods  in 
which  the  data  (d^}  is  smoothed  using  the  model 

n 

S (y)  = l a£nJ  4>.  (y)  , n < »,  (4.5) 

k=l  * * 


and  an  appropriate  statistical  fitting  strategy,  such  as  least  squares,  where 
(<t>k(y)}  belongs  to  an  explicitly  known  A-invertible  coordinate  system 
Consequently,  the  application  of  the  inversion  formula  is  accomplished 
analytically  to  yield  the  approximation 


2.  Solve  the  given  first  kind  Abel-type  equation  using  a stabilized 

computational  strategy  such  as  regularization  and,  if  necessary,  filter 
(smooth)  the  results.  Except  for  the  work  of  Baev  and  Glasko  (1976),  little 
appears  to  have  been  done  with  this  approach,  though  it  has  potential  when 
viewed  as  the  basis  for  a direct  filtering  strategy  in  the  sense  of  spectral 
differentiation. 


A third,  but  not  so  obvious  way,  is  to  perform  the  filtering  and 
stabilization  simultaneously.  This  is  the  basis  of  the  methods  proposed  by 
Andersscn  (1976)  and  studied  in  some  detail  by  Jakeman  (1975).  They  involve, 
through  the  use  of  product  integration  and  spectral  differentiation,  the 
direct  stabilization  of  the  inversion  formulas  (1.6a)  and  (1.6b)  so  that  they 
can  be  applied  to  the  given  discrete  observational  data  {d^ } . 

For  reasons  which  will  become  clearer  below,  we  shall  limit  attention 
to  an  examination  of  pseudo-analytic  methods  and  the  evaluation  of  stabilized 
inversion  formulas. 

Though  they  have  obvious  disadvantages  such  as 

die  (i)  The  smooth  representation  sn(y)  of  (4.5)  is  constructed  with 
a special  sense  (i.e.,  the  form  of  {4>g})  which  might  not 
be  representative  of  the  true  structure  in  the  data. 

die  (ii)  Since  (4.3)  is  weakly  improperly  posed,  the  use  of  the  smooth 
representation  sn(y)  which  does  not  truly  approximate  the 
structure  in  the  data  {d^}  could  yield  incorrect  results 
(e.g.,  see  Figure  2 in  Anderssen  (1973)). 

the  pseudo-analytic  methods  do  have  the  following  advantages 

ad  (i)  They  are  simple  to  implement  - once  a suitable  has 

been  chosen,  it  is  only  necessary  to  do  linear  parameter 
estimation  with  respect  to  an  appropriate  assumption  about 
the  structure  of  the  errors  in  the  data  {d^}. 

ad  (ii)  They  are  not  restricted  to  evenly  spaced  data. 

ad  (iii)  Assuming  that  sn(y)  of  (4.5)  is  a reasonable  approximation 

for  the  true  structure  in  the  data  and  the  parameter  estimation 
problem  is  numerically  stable,  then  the  totality  of  data  points 
required  for  a satisfactory  solution  is  potentially  not  too 
large. 


On  the  other  hand,  though  the  evaluation  of  stabilized  inversion  formulas  have 
obvious  advantages  which  more  or  less  correspond  to  the  negation  of  dis(i)  and 
die{ ii),  it  also  has  disadvantages,  directly  connected  with  the  use  of  spectral 
differentiation,  which  more  or  less  correspond  to  the  negation  of  ad(i)  , ac((ii) 
and  ad(iii) . 

Clearly,  the  key  point  is  that  they  potentially  define  independent 
rather  than  related  methods,  and  therefore  should  only  be  comoared  from  the  point 
of  view  of  identifying  their  compatibility.  For  example,  the  pseudo-analytic 
methods  can  be  viewed  as  simple  strategies  for  generating  first  order  approxima- 
tions to  the  required  solution  independently  of  the  signal  structure  in  the  data 
{d^},  but  not  more  accurate  approximations;  while  the  evaluation  of  the 
stabilized  inversion  formulas  defines  a sophisticated  method  for  computing  an 
accurate  approximation  to  the  required  solution  but  only  when  the  signal  structure 
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in  the  data  (d ^ } satisfies  conditions  compatible  with  such  a method  (i.e., 
with  the  use  of  spectral  differentiation). 

Thus,  the  sensible  strategy  would  appear  to  be  to  use  them  as 
compatible  methods  so  that  their  advantages  support  each  other  at  the  expense 
of  their  disadvantages.  In  particular,  for  given  data  {d^ } , generate  a 
first  approximation  to  the  required  solution  using  a pseudo-analytic  method 
in  such  a way  that  the  residual  data  has  a form  compatible  with  the  assumptions 
which  underlie  the  use  of  one  of  the  stabilized  inversion  formulas.  This  is 
the  basis  for  the  computational  strategy  proposed  by  Anderssen  (1976)  and 
studied  in  some  detail  in  Jakeman  (197S) 

More  explicit  details  about  the  algorithm  is  given  in  "Appendix  - The 
Algorithm". 


§5.  Linear  Functionals  and  the  Process  of  Inference 


r 


When  theoretical  and  numerical  investigations  of  special  classes  of 
equations  (such  as  Abel-type  integral  equations)  are  conducted  independently 
of  applications,  information  about  the  structure  of  the  inference  processes, 
which  use  the  solutions  generated  by  such  investigations  and  which  are  used 
to  make  conclusions  about  some  given  application,  are  ignored.  For  certain 
theoretical  investigations  (such  as  the  study  of  the  existence  and  smoothness 
of  solutions),  there  may  be  justification  for  ignoring  this  information.  But, 
as  we  shall  show  below  for  separable  Abel-type  equations  in  application,  this 
information  is  of  importance  for  some  aspects  of  theoretical  investigations 
(such  as  the  construction  of  inversion  formulas)  and  is  of  paramount  importance 
for  numerical  investigations,  as  it  may  allow  the  numerical  problem  to  be 
recast  in  a computationally  more  favourable  framework. 

A rationale  and  motivation  for  making  wherever  possible  such  a change 
of  framework  is:  Within  a more  favourable  computational  framework,  the 
probability  is  higher  that  inferences  consistent  with  the  information  in  given 
discrete  observational  data  will  be  generated. 

In  fact,  for  separable  Abel -type  integral  equations,  when  the  inference 
process  is  based  solely  on  the  use  of  linear  functionals  of  the  solution,  the 
mildly  improperly  posed  problem,  consisting  of  first  solving  the  relevant  Abel- 
type  equation  for  discrete  observational  data  followed  by  the  evaluation  of  the 
linear  functionals  on  the  resulting  approximate  solution,  can  be  replaced  by 
the  direct  evaluation  of  the  linear  functionals  on  the  given  discrete  observational 
data. 

To  illustrate,  we  consider  the  following  problem:  for  a given  consignment 
of  steel,  determine  whether  it  has  the  (impact)  strength  to  do  a given  job.  As 
the  result  of  extensive  experimentation,  metallurgists  claim  that  (Hyam  and 
Nutting  (1956)) 

(i)  Steel  consists  of  a conglomerate  of  ferrous  crystals  and  carbide 
(carbon-type)  particles,  which  to  a first  approximation  can  be 
regarded  as  spheres. 

(ii)  The  strength  of  steel  depends  on  its  hardness.  If  it  is  too  hard, 
it  is  too  brittle.  Tempering  is  used  to  make  it  softer  and  thereby 
increase  its  strength. 

(iii)  The  hardness  of  tempered  steels  depend  upon  the  ferrite  grain  sizes. 

The  general  rule  is  "the  larger  the  ferrite  grains,  the  softer  the 
steel". 

(iv)  The  growth  of  the  ferrite  grains  is  limited  by  the  carbide  particles. 

The  carbide  particles  form  along  grain  boundaries,  and  therefore  the 
size  of  the  ferrite  grains  is  heavily  constrained  by  the  size  distri- 
bution and  number  density  of  the  carbide  particles  (assuming  the 
centers  of  the  particles  are  distributed  something  like  a Poisson 
process) . 

(v)  During  tempering,  the  number  density  of  carbide  particles  decreases  but 
their  sizes  grow,  so  the  ferrite  grains  grow.  As  tempering  progresses, 
the  microstructure  and  mechanical  properties  continue  to  alter,  due  to 
the  development  of  approximately  spherical  particles  of  carbide 
(cemcntite)  which  gradually  grow  in  size  and  decrease  in  number. 
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The  mentioned  change  in  size  and  number  of  carbide  particles  can  only 
occur  by  the  solution  of  the  smaller  particles,  diffusion  of  carbon  through 
the  ferrite,  and  the  simultaneous  growth  of  the  larger  particles.  This 
indicates  that,  because  the  amount  of  ferrite,  carbon,  etc.  in  a given  consign- 
ment of  steel  being  tempered  is  fixed,  there  is  a trade-off  between  number 
density  p and  size  distribution  u(x)  of  carbide  particles  during  tempering 
at  the  expense  of  hardness  and  hence  at  the  gain  of  strength. 

Thus,  since  a knowledge  of  the  size  distribution  u(x)-  is  necessary 
for  estimating  the  number  density  p , the  problem  is  reduced  to  an  examination 
of  the  properties  of  the  size  distribution  u(x).  This  completes  the  first  step 
in  the  analysis  of  an  application  - pinpoint  the  properties  of  the  application 
which  relate  directly  to  the  problem  of  interest. 

The  second  step  is  to  derive  a mathematical  formalism  which  relates  u(x) 
to  observable  data.  As  a first  approximation,  we  shall  assume  that  the  carbide 
particles  are  spherical.  Since  it  is  not  possible  to  observe  the  carbide  spheres 
and  hence  u(x)  directly,  it  is  necessary  to  estimate  it  by  some  indirect  pro- 
cedure. It  is  because  of  this  necessity  to  observe  the  thing  one  requires  by 
some  indirect  procedure  that  improperly  posed  formulations  often  arise.  The  ways 
in  which  this  can  be  done  are  not  unique,  and  depend  heavily  on  the  sampling 
strategy  used.  The  one  which  we  shall  pursue  here  is  to  use  the  size  distribution 
of  circles  s(y)  (i.e.,  the  radius  or  diameter  distribution)  defined  on  random 
plane  sections  taken  through  the  steel.  The  sampling  strategy  is:  on  random 
plane  (polished)  sections  taken  through  a suitably  chosen  sample  of  the  steel, 
measure  the  observed  size  distribution  of  circles.  This  will  yield  the  data  {d^}. 

The  motivation  for  assuming  first  that  the  carbide  particle  «re  spherical 
(as  well  as  agreeing  approximately  with  observation)  and  for  taking  the  size 
distribution  of  circles  s(x),  as  the  means  by  which  we  estimate  u(x),  is  the 
existence  of  the  following  mathematical  formalism  (an  Abel-type  integral  equation) 
which  relates  u(x)  and  s (y) : 

a 

s(y)  = £ j u(x)  / (x2  - y2)'5  dx,  0 € y < x < a < <*>,  (5.1) 

y 

where  [0,a]  defines  the  support  of  u(x)  and 

a 

m = | xu(x)  dx  (5.2) 

0 

denotes  the  average  sphere  radius. 


Note  5.1.  The  actual  procedure,  for  choosing  a sample  of  steel  (from 
the  consignment),  and  then  taking  random  plane  sections  through  it,  is  non- 
trivial and  involves  a number  of  deep  (probabilistic)  questions.  A discussion 
of  some  of  the  difficulties  which  arise  along  with  proposals  for  using  volume 
weighted  and  area  weighted  sampling  can  be  found  in  the  references  of  §2.2. 

The  sampling  and  measurement  of  circles  will  yield  discrete  observations 
of  the  cumulative  distribution 


S(y)  = J s(z)  dz  , s(0)  = 0 , S(a)  = 1, 

i 0 

namely 
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di  = S(y.)  ♦ e.  . 0 < y,  < y2  < ....  < yn  . 

We  assume  that  the  measurement  process  is  sufficiently  accurate  so  that 
duplicate  observations  do  not  occur.  Once  n is  sufficiently  large  (so 
that  a sufficiently  detailed  amount  of  information  about  the  structure  of 
S(y)  has  been  accumulated),  the  third  step  can  be  implemented:  apply  a 
viable  computational  strategy  to  the  data  (d^ } to  generate  an  approximation 

{Gj}  = (u(Xj);  Xj  = (j  - 1)  Ax,  j = 1,  2,  ....  N,  Ax  = a/(N  - 1),  N < n) 

to  u(x)  on  an  even  grid  {xj}.  For  example,  apply  a two-stage  procedure  along 
the  lines  proposed  by  Anderssen  (1976),  implemented  by  Anderssen  and 
Jakeman  (1975b),  and  discussed  in  §4  above. 

The  fourth  step  involves  the  use  of  the  approximation  of  u(x) 

to  derive  inferences  about  the  strength  of  the  steel  in  the  given  consignment. 
But,  here  is  the  catch.  Though  the  metallurgist  argues  that  the  strength  of 
steel  is  governed  by  the  size  distribution  u(x)  of  the  carbide  particles,  he 
doesn't  necessarily  use  u(x)  directly  (see,  for  example,  Hilliard  (1968), 
p. 1373,  Underwood  (1970),  p.149,  Hyam  and  Nutting  (1956)}  when  making  a decision 
about  whether  a given  steel  has  the  strength  to  do  a given  job.  Instead,  he 
makes  inferences  on  the  basis  of  the  values  of  functionals  defined  on  u(x). 

The  functionals  of  u(x)  commonly  used  by  metallurgists  for  this  purpose  are 
(see,  for  example,  Hyam  and  Nutting  (1956)): 

a.  average  sphere  radius  m of  (5.2)  (which  we  assumed  above  to  be 
known) ; 


b.  average  particle  surface  area 

a 

Au  = 4tt  J x2u(x)dx  ; (5.3) 

0 

c.  average  particle  volume 

a 

vu  = f { x3u(x)dx  ; (S.4) 

0 


d.  number  of  particles  Ny  per  unit  volume  (assuming  their  centers  follow 
a Poisson  distribution) 


Ny  = NA  / (2m)  , (5.5) 

where  NA  denotes  the  number  of  circles  per  unit  area  of  planar  section 
of  the  sample  of  steel  on  which  the  data  (d^ } was  measured. 


Thus,  the  implementation  of  step  foul'  initially  involves  the  evaluation 
of  the  functionals  m,  Ay  and  V and  whence  Ny.  Invariably,  metallurgists 
estimate  m.  Ay  and  Vu  by  evaluating  the  integrals  (5.2),  (5.3)  and  (5.4) 
using  the  approximation  (u;)  - see,  for  example,  Hyam  and  Nutting  (1956), 

p.151,  column  2.  They  miss  the  point  that,  if  only  estimates  of  m.  Ay,  Vu  and 
Ny  are  required  for  inference  purposes,  then  the  necessity  to  explicitly  estimate 
u(x)  can  be  circumvented. 
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Because  inversion  formulas  for  separable  first  kind  Abel -type  integral 
equations  are  known  explicitly,  they  can  be  used  to  replace  linear  functionals 
defined  on  u(x),  namely 

F(0;z)  = | 8(x)u(x)dx, 

z 

by  linear  functionals  on  s(x),  namely 

a 

F(4>;z)  = | <Kz;x)  s(x)dx, 
z 

In  fact,  for  a given  9(x),  the  4>(z;x)  corresponding  to  the  Abel  integral 
equation  (4.3),  with  kj(y)  = y/m  and  k2(x)  = 1,  is  given  by 

x 

*(z;x)  = ^ (6(z) (x2  - z2)'h  * J ^ (x2  - y*)~hdy},  (5.6) 

z 

lA  1 

if  we  assume  that  0(x)  is  such  that  ^-(x2  - y2)‘^  is  integrable,  and  0(y)  f y. 

When  0(y)  = y,  it  is  necessary  to  employ  the  type  of  independent  argument  used 
by  Watson  (1971),  since  the  above  formula  (5.6)  for  <t>(x)  involves  m.  The 
existence  of  the  general  correspondence  (5.6)  between  0(x)  and  4> (x)  appear  to 
have  been  overlooked  in  practice,  except  for  the  special  cases  examined  by 
Meisner  (1967),  Watson  (1971)  and  Jakeman  (1975;  eqn(4.3.1)). 

Thus,  the  problem  of  estimating  the  functionals  (5.2),  (5.3),  (5.4)  and 
(5.5)  reduces  to  evaluating  the  functionals 

m = \ ( | ^ dx  )'*  (5.7) 

0 

a 

Au  = 16m  | xs(x)dx  (5.8) 

0 

a 

Vu  = 2trm  | x2s(x)dx  (5.9) 

0 

Ny  = NA  / (2m)  (5.10) 

which  are  properly  posed  formulations.  Though  standard  quadrature  rules,  such 
as  the  rectangular,  trapezoidal  and  Simpson,  can  be  applied  to  evaluate  the 
above  functionals  for  the  data  {di> , it  doesn't  follow  that  they  are  necessarily 
the  most  appropriate.  Not  only  is  it  necessary  to  choose  a quadrature  which 
performs  well  numerically;  it  is  also  necessary,  because  of  the  observational 
nature  of  the  data  {dA ) {see  eqn(4.4)},  to  choose  quadrature  rules  which  yield 
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estimators  of  the  above  functionals  which  have  appropriate  statistical  properties. 
This  problem  has  been  examined  in  some  detail  by  Anderssen  and  Jakeman  (1975a) 
and  Jakeman  and  Scheaffer  (1977).  They  conclude  that,  especially  for  functionals 
involving  integrals  of  the  form 


1 


_ s(x) 

(x2-z2?* 


dx. 


0 <■  z < a. 


product  integration  procedures  should  be  used. 


From  a general  point  of  view,  the  above  discussion  illustrates  that  the 
strategy  for  problem  solving  is  not  simply 

(i)  State  (define)  the  problem. 

(ii)  Identify  the  appropriate  (experimental,  observational  and/or 
theoretical)  facts  about  the  problem  and  its  context  from  which  one 
derives  structural  assumptions  for  the  construction  of  a framework 
into  which  the  problem  can  be  embedded  rigorously  (mathematically, 
scientifically,  logically  or  unambiguously). 

(iii)  Define  a mathematical  formalism  for  the  problem,  and,  when 
necessary,  an  appropriate  sampling  strategy  for  the  data. 

(iv)  Through  the  use  of  an  appropriate  rationale  (consistent  with 
the  mathematical  formalism),  derive  a viable  computational  formalism 
for  the  mathematical  formalism. 

(v)  Construct  an  algorithmic  implementation  for  the  computational 
formalism,  and  evaluate  a specific  approximation  to  the  required 
solution  of  the  mathematical  formalism. 

(vi)  Use  the  resulting  approximation  for  inference  purposes  about 
the  problem  context. 

In  fact,  the  strategy  should  be  (assuming  backtracking  takes  place 
when  and  where  necessary) 

a.  (i) 

b.  (ii) 

c.  (iii) 

d.  Identify  the  functionals,  defined  on  the  solution  of  the 

mathematical  formalism,  from  which  inferences  about  the  problem  will 
be  derived. 

e.  Through  the  use  of  an  appropriate  rationale  (consistent  with 
the  mathematical  formalism  and  the  identified  functionals),  derive  a 
viable  computational  formalism  for  the  mathematical  formalism. 

f.  Construct  an  algorithmic  implementation  for  the  computational 

formalism,  and  evaluate  a specific  approximation  to  the  required 
functionals. 

g.  Use  the  resulting  estimates  of  the  functionals  for  inference 
purposes. 
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APPENDIX  - TABLE  1 (continued) 


lPPENDIX  - TABLE  1 (continued) 
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Table  of  K(a,x),  Corresponding  Form  of  (3.1),  and  Associated  Site  Distribution  Estimation  Problem 
Data  Site  Distribution  Estimation  Problem  K(a,x)  Corresponding  Form  of  Equation  (3. 
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An  algorithm  has  been  implemented  and  copies  of  the  FORTRAN  program  are 
available  from  either  R.  S.  Anderssen  (Computer  Centre)  or  A.  J.  Jakeman  (CRES) , 
Australian  National  University,  P.  0.  Box  4,  Canberra  ACT2600,  Australia.  At 
the  montent,  it  is  only  implemented  for  the  random  spheres  model  (see  Appendix  - 
Table  1,  application  lb).  A general  package  for  Abel-type  integral  equations 
along  with  test  data  is  in  preparation.  It  will  include  algorithms  for  the  direct 
evaluation  of  linear  functionals  defined  on  the  data. 
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